MODULE VEXCSSAD_MOD
CONTAINS
SUBROUTINE VEXCSSAD (KIDIA,KFDIA,KLON,PTMST,PRVDIFTS, &
 & LDEKPERTS, LDREGBUOF, &
 & PUMLEV5,PVMLEV5  ,PTMLEV5   ,PQMLEV5, &
 & PAPHMS5,PGEOMLEV5,PCPTGZLEV5,PCPTS5 , &
 & PQSAM5 ,PZ0MM5   ,PZ0HM5    ,PZ0QM5 , PBUOM5 , &
 & PUCURR5,PVCURR5  , &
 & YDCST  ,YDEXC    , &
 & PCFM5  ,PCFH5    ,PCFQ5     , &
 & PUMLEV ,PVMLEV   ,PTMLEV    ,PQMLEV , &
 & PAPHMS ,PGEOMLEV ,PCPTGZLEV ,PCPTS  , &
 & PQSAM  ,PZ0MM    ,PZ0HM     ,PZ0QM  , PBUOM  , &
 & PCFM   ,PCFH     ,PCFQ )

USE PARKIND1  , ONLY : JPIM, JPRB
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_THF   , ONLY : RVTMP2
USE YOS_EXCS  , ONLY : RLPDD, RLPCC, RLPBB
USE YOS_CST   , ONLY : TCST
USE YOS_EXC   , ONLY : TEXC

! (C) Copyright 1995- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ------------------------------------------------------------------

!**   *VEXCSSAD* - DETERMINES THE EXCHANGE COEFFICIENTS BETWEEN THE
!                   SURFACE AND THE LOWEST MODEL LEVEL WITH HELP OF
!                   STABILITY AS FUNCTION OF RICHARDSON NUMBER.
!                   (Adjoint)

!     J.F. MAHFOUF                E.C.M.W.F.    02/10/95
!     Modified   P. Viterbo ECMWF 12/05/2005 Externalize SURF
!                                   (based on vdfexcstl)
!                M. Janiskova  29/03/2006  Code optimization and cleaning
!                M. Janiskova  08/08/2007  Changed regularization
!                M. Janiskova  Aug 2011    Adjusted regularization for
!                                          momentum exch. coefficient
!                P. Lopez      July 2025   Added ocean currents (trajectory only)
!                P. Lopez      July 2025   Added optional (LDREGBUOF) extra regularization 
!                                          when surface buoyancy flux is very small.

!     PURPOSE
!     -------

!     DETERMINE EXCHANGE COEFFICIENTS BETWEEN THE SURFACE AND THE
!     LOWEST MODEL LEVEL

!     INTERFACE
!     ---------

!     *VEXCSTL* IS CALLED BY *SURFEXCDRIVERTL*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLON*         NUMBER OF GRID POINTS PER PACKET

!     INPUT PARAMETERS (LOGICAL)

!     LDEKPERTS      TRUE IF PERTURBATION OF EXCHANGE COEEFICIENTS
!     LDREGBUOF      TRUE FOR EXTRA REGULARIZATION WHEN SURFACE BUOYANCY FLUX IS VERY SMALL

!     INPUT PARAMETERS (REAL):

!     *PTMST*        TIME STEP
!     *PRVDIFTS*     Semi-implicit factor for vertical diffusion discretization

!  Trajectory  Perturbation  Description                               Unit
!  PUMLEV5     PUMLEV        X-VELOCITY COMPONENT AT T-1, 
!                            lowest model level                        m/s
!  PVMLEV5     PVMLEV        Y-VELOCITY COMPONENT AT T-1, 
!                            lowest model level                        m/s
!  PTMLEV5     PTMLEV        TEMPERATURE AT T-1, lowest model level    K
!  PQMLEV5     PQMLEV        SPECIFIC HUMUDITY AT T-1, 
!                            lowest model level                        kg/kg
!  PAPHMS5     PAPHMS        PRESSURE AT T-1 , surface                 Pa
!  PGEOMLEV5   PGEOMLEV      GEOPOTENTIAL AT T-1, lowest model level   m2/s2
!  PCPTGZLEV5  PCPTGZLEV     DRY STATIC ENERGY AT LOWEST MODEL LEVEL   J/kg
!  PCPTS5      PCPTS         DRY STATIC ENERGY AT THE SURFACE          J/kg
!  PQSAM5      PQSAM         SPECIFIC HUMIDITY AT THE SURFACE          kg/kg
!  PZ0MM5      PZ0MM         AERODYNAMIC ROUGHNESS LENGTH              m
!  PZ0HM5      PZ0HM         ROUGHNESS LENGTH FOR TEMPERATURE          m
!  PZ0QM5      PZ0QM         ROUGHNESS LENGTH FOR MOISTURE             m
!  PBUOM5      PBUOM         BUOYANCY FLUX AT THE SURFACE              ?
!  PUCURR5     -----         OCEAN CURRENT U-COMPONENT                 m/s
!  PVCURR5     -----         OCEAN CURRENT V-COMPONENT                 m/s

!     OUTPUT PARAMETERS (REAL):

!  Trajectory  Perturbation  Description                               Unit
!  PCFM5       PCFM          PROP. TO EXCH. COEFF. FOR MOMENTUM
!                            (C-STAR IN DOC.) AT LOWEST MODEL LEVEL    ?
!  PCFH5       PCFH          PROP. TO EXCH. COEFF. FOR HEAT
!                            (C-STAR IN DOC.) AT LOWEST MODEL LEVEL    ?
!  PCFQ5       PCFQ          PROP. TO EXCH. COEFF. FOR MOISTURE
!                            (C-STAR IN DOC.) AT LOWEST MODEL LEVEL    ?

!     METHOD
!     ------

!     STABILITY FUNCTIONS ARE EXPLICIT FUNCTIONS OF THE BULK RICHARDSON
!     NUMBER ACCORDING TO LOUIS ET AL. (1982) AND GENERALIZED TO THE
!     SITUATION OF DIFFERENT ROUGHNESS LENGTHS FOR HEAT AND MOMENTUM

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
LOGICAL, INTENT(IN)              :: LDEKPERTS
LOGICAL, INTENT(IN)              :: LDREGBUOF
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMST 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRVDIFTS
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGEOMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTGZLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTS5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQSAM5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0MM5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0HM5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0QM5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PBUOM5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUCURR5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVCURR5(:) 
TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TEXC)        ,INTENT(IN)    :: YDEXC
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFM5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFH5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFQ5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PUMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PVMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PTMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PAPHMS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PGEOMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCPTGZLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCPTS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQSAM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PZ0MM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PZ0HM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PZ0QM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBUOM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFH(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFQ(:) 

!*            LOCAL STORAGE
!             ----- -------

REAL(KIND=JPRB) ::  Z1DZ0M5(KLON),Z1DZ0H5(KLON), Z1DZ0Q5(KLON) ,ZXLNM5(KLON) ,&
 & ZXLNH5(KLON) ,ZXLNQ5(KLON) , ZDU25(KLON)   ,ZNLEV5(KLON)   
REAL(KIND=JPRB) :: ZRICLS5(KLON), ZQSAM5(KLON)

REAL(KIND=JPRB) ::  Z1DZ0M(KLON)  ,Z1DZ0H(KLON),  Z1DZ0Q(KLON)  ,ZXLNM(KLON) ,&
 & ZXLNH(KLON)   ,ZXLNQ(KLON) ,  ZDU2(KLON)    ,ZNLEV(KLON)
REAL(KIND=JPRB) :: ZRICLS(KLON), ZQSAM(KLON)

INTEGER(KIND=JPIM) :: JL

REAL(KIND=JPRB) :: Z1S, Z1S5, Z2B, Z2S, Z2S5, Z3B, Z5S, Z5S5,&
 & Z6S, Z6S5, Z7S, Z7S5, ZAUX1, ZAUX15, ZAUX2, &
 & ZAUX25, ZC, ZCDNH, ZCDNH5, ZCDNM, ZCDNM5, &
 & ZCDNQ, ZCDNQ5, ZCH, ZCH5, ZCM, ZCM5, ZCON1, &
 & ZCON2, ZCONS12, ZCQ, ZCQ5, ZD, ZDRORO, ZDRORO5, &
 & ZIPBL, ZRIMIN, ZTPFAC1, ZUABS, ZUABS5, ZWST2, &
 & ZWST25, ZRGI, ZRKAP3B
REAL(KIND=JPRB) :: ZDIV15, ZDIV25, ZDIV35, ZDIV45, ZDIV55, ZDIV65
REAL(KIND=JPRB) :: ZDIV75, ZDIV85, ZDIV95, ZDIV105, ZDIV115
REAL(KIND=JPRB) :: ZCFQ5, ZCFH5, ZCFM5
REAL(KIND=JPRB) :: ZAA, ZCC, ZAA1, ZCC1, ZDD, ZDD1
REAL(KIND=JPRB) :: ZAA2, ZCC2, ZDD2
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

!*       1.   INITIALIZE CONSTANTS
!             ---------- ----------

IF (LHOOK) CALL DR_HOOK('VEXCSSAD_MOD:VEXCSSAD',0,ZHOOK_HANDLE)
ASSOCIATE(RD=>YDCST%RD, RETV=>YDCST%RETV, RG=>YDCST%RG, &
 & REPDU2=>YDEXC%REPDU2, RKAP=>YDEXC%RKAP, RPARZI=>YDEXC%RPARZI)
ZTPFAC1 = PRVDIFTS
ZCONS12 = (ZTPFAC1*PTMST*RG)/RD
ZCON1 = RVTMP2-RETV
ZCON2 = 2.0_JPRB/3._JPRB
ZRIMIN = 1.E-3_JPRB
ZRGI = 1.0_JPRB/RG

!     CONSTANTS FOR REGULARIZATION

ZAA=5.0_JPRB
ZCC=0.0416277_JPRB
ZAA1=20.0_JPRB
ZCC1=0.0588705_JPRB
ZAA2=10.0_JPRB
ZCC2=0.0416277_JPRB

!     CONSTANTS FOR THE STABILITY FUNCTIONS

Z2B = 2.0_JPRB*RLPBB
ZD = RLPDD
ZC = RLPCC
Z3B = 3._JPRB*RLPBB
ZRKAP3B = RKAP*RKAP*Z3B

!     PBL HEIGHT FOR W* - EFFECT

ZIPBL = RPARZI

!------------------------------------------------------------------------

!*                ------- TRAJECTORY COMPUTATIONS ---------

!------------------------------------------------------------------------

!*       2.   COMPUTATION OF BASIC QUANTITIES
!             ----------- -- ----- ----------

DO JL = KIDIA, KFDIA

!             (W*)**2, WIND SHEAR, RICHARDSON NUMBER,

  IF (PBUOM5(JL) <= 0.0_JPRB) THEN
    ZWST25 = 0.0_JPRB
  ELSE
    ZWST25 = (ZIPBL*PBUOM5(JL))** ZCON2
  ENDIF
  Z1S5 = (PUMLEV5(JL)-PUCURR5(JL))**2 + (PVMLEV5(JL)-PVCURR5(JL))**2 + ZWST25
  IF (REPDU2 >= Z1S5) THEN
    ZDU25(JL) = REPDU2
  ELSE
    ZDU25(JL) = Z1S5
  ENDIF 
  ZDIV15 = 1.0_JPRB/(PCPTGZLEV5(JL)+PCPTS5(JL)-PGEOMLEV5(JL))
  IF (PQSAM5(JL) < 0.1_JPRB) THEN
    ZQSAM5(JL) = PQSAM5(JL)
  ELSE
    ZQSAM5(JL) = 0.1_JPRB
  ENDIF
  ZDRORO5 = (2.0_JPRB*(PCPTGZLEV5(JL)-PCPTS5(JL)))*ZDIV15 &
   & -ZCON1*(PQMLEV5(JL)-ZQSAM5(JL))  
  ZDIV25 = 1.0_JPRB/ZDU25(JL)
  ZRICLS5(JL) =  PGEOMLEV5(JL)*ZDRORO5*ZDIV25 

!             COMMON FACTORS IN NEUTRAL FORMULAE AND
!             DRAG COEFFICIENTS.

  ZNLEV5(JL) = PGEOMLEV5(JL)*ZRGI+PZ0MM5(JL)
  ZDIV35 = 1.0_JPRB/PZ0MM5(JL)
  Z1DZ0M5(JL) = ZNLEV5(JL)*ZDIV35
  ZDIV45 = 1.0_JPRB/PZ0HM5(JL)
  Z1DZ0H5(JL) = ZNLEV5(JL)*ZDIV45
  ZDIV55 = 1.0_JPRB/PZ0QM5(JL)
  Z1DZ0Q5(JL) = ZNLEV5(JL)*ZDIV55
  ZXLNM5(JL) = LOG(Z1DZ0M5(JL))
  ZXLNH5(JL) = LOG(Z1DZ0H5(JL))
  ZXLNQ5(JL) = LOG(Z1DZ0Q5(JL))

  ZCDNH5 = 1.0_JPRB/(ZXLNM5(JL)*ZXLNH5(JL))
  ZCDNQ5 = 1.0_JPRB/(ZXLNM5(JL)*ZXLNQ5(JL))
  ZCDNM5 = 1.0_JPRB/(ZXLNM5(JL)*ZXLNM5(JL))
  ZCM5 = ZC*SQRT(Z1DZ0M5(JL))
  ZCH5 = ZC*SQRT(Z1DZ0H5(JL))
  ZCQ5 = ZC*SQRT(Z1DZ0Q5(JL))

!           STABILITY FUNCTIONS 

  IF (ZRICLS5(JL) > 0) THEN
    Z2S5 = SQRT(1.0_JPRB+ZD*ZRICLS5(JL))
    ZDIV65 = 1.0_JPRB/(Z2S5+Z2B*ZRICLS5(JL))
    PCFM5(JL) = ZCDNM5*Z2S5*ZDIV65
    ZDIV75 = 1.0_JPRB/(1.0_JPRB+Z3B*ZRICLS5(JL)*Z2S5)
    PCFH5(JL) = ZCDNH5*ZDIV75
    PCFQ5(JL) = ZCDNQ5*ZDIV75
  ELSE

!  Threshold for minimum Richardson number

    IF (-ZRICLS5(JL) > ZRIMIN) THEN
      Z2S5 = SQRT(-ZRICLS5(JL))
    ELSE
      Z2S5 = SQRT(ZRIMIN)
    ENDIF
    Z5S5 = ZRKAP3B*ZCDNM5*ZCM5
    ZDIV85 = 1.0_JPRB/(1.0_JPRB+Z5S5*Z2S5) 
    PCFM5(JL) = ZCDNM5*(1.0_JPRB-(Z2B*ZRICLS5(JL))*ZDIV85)
    Z6S5 = ZRKAP3B*ZCDNH5*ZCH5
    ZDIV95 = 1.0_JPRB/(1.0_JPRB+Z6S5*Z2S5)
    PCFH5(JL) = ZCDNH5*(1.0_JPRB-(Z3B*ZRICLS5(JL))*ZDIV95)
    Z7S5 = ZRKAP3B*ZCDNQ5*ZCQ5
    ZDIV105 = 1.0_JPRB/(1.0_JPRB+Z7S5*Z2S5)
    PCFQ5(JL) = ZCDNQ5*(1.0_JPRB-(Z3B*ZRICLS5(JL))*ZDIV105)
  ENDIF
  ZDIV115 = 1.0_JPRB/(PTMLEV5(JL)*(1.0_JPRB+RETV*PQMLEV5(JL)))
  ZAUX15 = (ZCONS12*PAPHMS5(JL))*ZDIV115
  ZUABS5 = SQRT(ZDU25(JL))
  ZAUX25 = ZAUX15*ZUABS5*RKAP**2
  ZCFM5 = PCFM5(JL)
  PCFM5(JL) = ZAUX25*PCFM5(JL)
  ZCFH5 = PCFH5(JL)
  PCFH5(JL) = ZAUX25*PCFH5(JL)
  ZCFQ5 = PCFQ5(JL)
  PCFQ5(JL) = ZAUX25*PCFQ5(JL)

!------------------------------------------------------------------------

!*                ------- ADJOINT COMPUTATIONS ---------

!------------------------------------------------------------------------

!*                          INITIALISATIONS
!                           ---------------

  ZXLNQ  (JL) = 0.0_JPRB
  ZXLNH  (JL) = 0.0_JPRB
  ZXLNM  (JL) = 0.0_JPRB
  ZDU2   (JL) = 0.0_JPRB
  Z1DZ0Q (JL) = 0.0_JPRB
  Z1DZ0H (JL) = 0.0_JPRB
  Z1DZ0M (JL) = 0.0_JPRB
  ZNLEV  (JL) = 0.0_JPRB
  ZRICLS (JL) = 0.0_JPRB
  ZQSAM  (JL) = 0.0_JPRB

  ZAUX2  = 0.0_JPRB
  ZUABS  = 0.0_JPRB
  ZAUX1  = 0.0_JPRB
  ZCDNQ  = 0.0_JPRB
  Z7S    = 0.0_JPRB
  Z2S    = 0.0_JPRB
  ZCQ    = 0.0_JPRB
  ZCDNH  = 0.0_JPRB
  Z6S    = 0.0_JPRB
  ZCH    = 0.0_JPRB
  ZCDNM  = 0.0_JPRB
  Z5S    = 0.0_JPRB
  ZCM    = 0.0_JPRB
  ZWST2  = 0.0_JPRB
  Z1S    = 0.0_JPRB
  ZDRORO = 0.0_JPRB

  IF (.NOT.LDEKPERTS) THEN
    PCFM(JL) = 0.0_JPRB
    PCFH(JL) = 0.0_JPRB
    PCFQ(JL) = 0.0_JPRB
  ELSE
    IF (LDREGBUOF .AND. ABS(PBUOM5(JL)) < 2.E-6_JPRB) THEN
      PCFQ(JL) = PCFQ(JL) * 0.01_JPRB
      PCFH(JL) = PCFH(JL) * 0.01_JPRB
      PCFM(JL) = PCFM(JL) * 0.01_JPRB
    ENDIF

    ZDD  = MAX(ZAA*EXP(-(ZCC*ZRICLS5(JL))**2),1.0_JPRB)
    ZDD  = 1.0_JPRB/ZDD
    ZDD1 = MAX(ZAA1*EXP(-(ZCC1*ZRICLS5(JL))**2),1.0_JPRB)
    ZDD1 = 1.0_JPRB/ZDD1 
    ZDD2 = MAX(ZAA2*EXP(-(ZCC2*ZRICLS5(JL))**2),1.0_JPRB)
    ZDD2 = 1.0_JPRB/ZDD2

    IF (PTMLEV5(JL) > 263.15_JPRB) THEN
      PCFQ(JL) = PCFQ(JL)*ZDD
    ELSE
      PCFQ(JL) = PCFQ(JL)*ZDD1
    ENDIF
    IF (PTMLEV5(JL) > 270.65_JPRB .AND. ZUABS5 > 6.0_JPRB) THEN
      IF (PZ0HM5(JL) > 0.3E-04_JPRB .AND. PCFH5(JL) < 60.0_JPRB) THEN
        PCFH (JL) = 0.05_JPRB*PCFH (JL)
      ENDIF
      IF (PZ0HM5(JL) > 0.7E-03_JPRB) THEN
        PCFH (JL) = 0.05_JPRB*PCFH (JL)
      ENDIF
    ELSE
      PCFH (JL) = PCFH (JL)*ZDD1
    ENDIF
    IF (PTMLEV5(JL) < 270.65_JPRB .AND. ZUABS5 < 6.0_JPRB) THEN
      PCFM(JL) = PCFM(JL)*ZDD
    ELSEIF (ZUABS5 < 1.0_JPRB) THEN
      PCFM(JL) = PCFM(JL)*ZDD2
    ENDIF
  ENDIF 

!             STABILITY FUNCTIONS

  ZAUX2 =  ZAUX2 + PCFQ(JL)*ZCFQ5
  PCFQ(JL) = ZAUX25*PCFQ(JL)
  ZAUX2 =  ZAUX2 + PCFH(JL)*ZCFH5
  PCFH(JL) = ZAUX25*PCFH(JL)
  ZAUX2 =  ZAUX2 + PCFM(JL)*ZCFM5
  PCFM(JL) = ZAUX25*PCFM(JL)

  ZUABS = ZUABS + ZAUX2*ZAUX15*RKAP**2
  ZAUX1 = ZAUX1 + ZAUX2*ZUABS5*RKAP**2

  ZDU2 (JL) = ZDU2(JL)+0.5_JPRB*ZUABS/ZUABS5

  PAPHMS(JL) = PAPHMS(JL)+ZCONS12*ZAUX1*ZDIV115
  PQMLEV(JL) = PQMLEV(JL)-ZCONS12*ZAUX1*PAPHMS5(JL)*PTMLEV5(JL) &
   & *RETV*ZDIV115*ZDIV115
  PTMLEV(JL) = PTMLEV(JL)-ZCONS12*ZAUX1*PAPHMS5(JL) &
   & *(1.0_JPRB+RETV*PQMLEV5(JL))*ZDIV115*ZDIV115

  IF (ZRICLS5(JL) > 0) THEN
    ZCDNQ = ZCDNQ+ZDIV75*PCFQ(JL)
    Z2S = Z2S-Z3B*ZCDNQ5*ZDIV75*ZDIV75*ZRICLS5(JL)*PCFQ(JL)
    ZRICLS(JL) = ZRICLS(JL)-Z3B*ZCDNQ5*ZDIV75*ZDIV75*Z2S5*PCFQ(JL)
    PCFQ(JL) = 0.0_JPRB

    ZCDNH = ZCDNH+ZDIV75*PCFH(JL)
    Z2S = Z2S-Z3B*ZCDNH5*ZDIV75*ZDIV75*ZRICLS5(JL)*PCFH(JL)
    ZRICLS(JL) = ZRICLS(JL)-Z3B*ZCDNH5*ZDIV75*ZDIV75*Z2S5*PCFH(JL)
    PCFH(JL) = 0.0_JPRB

    ZCDNM = ZCDNM+Z2S5*ZDIV65*PCFM(JL)
    Z2S = Z2S+Z2B*ZCDNM5*ZDIV65*ZDIV65*ZRICLS5(JL)*PCFM(JL)
    ZRICLS(JL) = ZRICLS(JL)-Z2B*ZCDNM5*ZDIV65*ZDIV65*Z2S5*PCFM(JL)
    PCFM(JL) = 0.0_JPRB

    ZRICLS(JL) = ZRICLS(JL)+0.5_JPRB*ZD*Z2S/Z2S5
  ELSE
    ZCDNQ = ZCDNQ+(1.0_JPRB-(Z3B*ZRICLS5(JL))*ZDIV105)*PCFQ(JL)
    ZRICLS(JL) = ZRICLS(JL)-ZCDNQ5*Z3B*ZDIV105*PCFQ(JL)
    Z7S = Z7S+ZCDNQ5*Z3B*Z2S5*ZRICLS5(JL)*ZDIV105*ZDIV105*PCFQ(JL)
    Z2S = Z2S+ZCDNQ5*Z3B*Z7S5*ZRICLS5(JL)*ZDIV105*ZDIV105*PCFQ(JL)
    PCFQ(JL) = 0.0_JPRB
    ZCDNQ = ZCDNQ+ZRKAP3B*ZCQ5*Z7S
    ZCQ = ZCQ+ZRKAP3B*ZCDNQ5*Z7S

    ZCDNH = ZCDNH+(1.0_JPRB-(Z3B*ZRICLS5(JL))*ZDIV95)*PCFH(JL)
    ZRICLS(JL) =ZRICLS(JL)-ZCDNH5*Z3B*ZDIV95*PCFH(JL)
    Z6S = Z6S+ZCDNH5*Z3B*Z2S5*ZRICLS5(JL)*ZDIV95*ZDIV95*PCFH(JL)
    Z2S = Z2S+ZCDNH5*Z3B*Z6S5*ZRICLS5(JL)*ZDIV95*ZDIV95*PCFH(JL)
    PCFH(JL) = 0.0_JPRB
    ZCDNH = ZCDNH+ZRKAP3B*ZCH5*Z6S
    ZCH = ZCH+ZRKAP3B*ZCDNH5*Z6S

    ZCDNM = ZCDNM+(1.0_JPRB-(Z2B*ZRICLS5(JL))*ZDIV85)*PCFM(JL)
    ZRICLS(JL) = ZRICLS(JL)-ZCDNM5*Z2B*ZDIV85*PCFM(JL)
    Z5S = Z5S+ZCDNM5*Z2B*Z2S5*ZRICLS5(JL)*ZDIV85*ZDIV85*PCFM(JL)
    Z2S = Z2S+ZCDNM5*Z2B*Z5S5*ZRICLS5(JL)*ZDIV85*ZDIV85*PCFM(JL)
    PCFM(JL) = 0.0_JPRB
    ZCDNM = ZCDNM+ZRKAP3B*ZCM5*Z5S
    ZCM = ZCM+ZRKAP3B*ZCDNM5*Z5S

!  Threshold for minimum Richardson number

    IF (-ZRICLS5(JL) > ZRIMIN) THEN
      ZRICLS(JL) = ZRICLS(JL)-0.5_JPRB*Z2S/Z2S5
    ENDIF
  ENDIF

!             COMMON FACTORS IN NEUTRAL FORMULAE AND
!             DRAG COEFFICIENTS.

  Z1DZ0Q(JL) = Z1DZ0Q(JL)+0.5_JPRB*ZC*ZC*ZCQ/ZCQ5
  Z1DZ0H(JL) = Z1DZ0H(JL)+0.5_JPRB*ZC*ZC*ZCH/ZCH5
  Z1DZ0M(JL) = Z1DZ0M(JL)+0.5_JPRB*ZC*ZC*ZCM/ZCM5

  ZXLNM(JL) = ZXLNM(JL)-2*ZXLNM5(JL)*ZCDNM5*ZCDNM5*ZCDNM
  ZXLNQ(JL) = ZXLNQ(JL)-ZXLNM5(JL)*ZCDNQ5*ZCDNQ5*ZCDNQ
  ZXLNM(JL) = ZXLNM(JL)-ZXLNQ5(JL)*ZCDNQ5*ZCDNQ5*ZCDNQ
  ZXLNH(JL) = ZXLNH(JL)-ZXLNM5(JL)*ZCDNH5*ZCDNH5*ZCDNH
  ZXLNM(JL) = ZXLNM(JL)-ZXLNH5(JL)*ZCDNH5*ZCDNH5*ZCDNH

  Z1DZ0Q(JL) = Z1DZ0Q(JL)+ZXLNQ(JL)/Z1DZ0Q5(JL)
  ZXLNQ(JL) = 0.0_JPRB
  Z1DZ0H(JL) = Z1DZ0H(JL)+ZXLNH(JL)/Z1DZ0H5(JL)
  ZXLNH(JL) = 0.0_JPRB
  Z1DZ0M(JL) = Z1DZ0M(JL)+ZXLNM(JL)/Z1DZ0M5(JL)
  ZXLNM(JL) = 0.0_JPRB
  ZNLEV(JL) = ZNLEV(JL)+PZ0QM5(JL)*ZDIV55*ZDIV55*Z1DZ0Q(JL)
  PZ0QM(JL) = PZ0QM(JL)-ZNLEV5(JL)*ZDIV55*ZDIV55*Z1DZ0Q(JL)
  Z1DZ0Q(JL) = 0.0_JPRB
  ZNLEV(JL) = ZNLEV(JL)+PZ0HM5(JL)*ZDIV45*ZDIV45*Z1DZ0H(JL)
  PZ0HM(JL) = PZ0HM(JL)-ZNLEV5(JL)*ZDIV45*ZDIV45*Z1DZ0H(JL)
  Z1DZ0H(JL) = 0.0_JPRB
  ZNLEV(JL) = ZNLEV(JL)+PZ0MM5(JL)*ZDIV35*ZDIV35*Z1DZ0M(JL)
  PZ0MM(JL) = PZ0MM(JL)-ZNLEV5(JL)*ZDIV35*ZDIV35*Z1DZ0M(JL)
  Z1DZ0M(JL) = 0.0_JPRB
  PGEOMLEV(JL) = PGEOMLEV(JL)+ZRGI*ZNLEV(JL)
  PZ0MM(JL) = PZ0MM(JL)+ZNLEV(JL)
  ZNLEV(JL) = 0.0_JPRB

!             (W*)**2, WIND SHEAR, RICHARDSON NUMBER,

  PGEOMLEV(JL) = PGEOMLEV(JL)+ZDRORO5*ZDIV25*ZRICLS(JL)
  ZDRORO = ZDRORO+PGEOMLEV5(JL)*ZDIV25*ZRICLS(JL)
  ZDU2(JL) = ZDU2(JL)-PGEOMLEV5(JL)*ZDRORO5*ZDIV25*ZDIV25*ZRICLS(JL)
  ZRICLS(JL) = 0.0_JPRB
  PGEOMLEV(JL) = PGEOMLEV(JL)+(2.0_JPRB*(PCPTGZLEV5(JL)-PCPTS5(JL))) &
   & *ZDIV15*ZDIV15*ZDRORO
  PCPTGZLEV(JL) = PCPTGZLEV(JL)+2.0_JPRB &
   & *(2.0_JPRB*PCPTS5(JL)-PGEOMLEV5(JL))*ZDIV15*ZDIV15*ZDRORO
  PCPTS(JL) = PCPTS(JL)-2.0_JPRB &
   & *(2.0_JPRB*PCPTGZLEV5(JL)-PGEOMLEV5(JL))*ZDIV15*ZDIV15*ZDRORO
  PQMLEV(JL) = PQMLEV(JL)-ZCON1*ZDRORO
  ZQSAM(JL) = ZQSAM(JL)+ZCON1*ZDRORO
  IF (PQSAM5(JL) < 0.1_JPRB) THEN
    PQSAM(JL) = PQSAM(JL)+ZQSAM(JL)
  ENDIF
  ZQSAM(JL) = 0.0_JPRB
  IF (REPDU2 >= Z1S5) THEN
    ZDU2 (JL) = 0.0_JPRB
  ELSE
    Z1S = Z1S+ZDU2(JL)
    ZDU2 (JL) = 0.0_JPRB
  ENDIF
  PUMLEV(JL) = PUMLEV(JL) + 2.0_JPRB*(PUMLEV5(JL)-PUCURR5(JL))*Z1S
  PVMLEV(JL) = PVMLEV(JL) + 2.0_JPRB*(PVMLEV5(JL)-PVCURR5(JL))*Z1S
  ZWST2 = ZWST2+Z1S
  IF (PBUOM5(JL) <= 0.0_JPRB) THEN
    ZWST2  = 0.0_JPRB
  ELSE
    PBUOM(JL) = PBUOM(JL)+ZCON2*ZIPBL*ZWST2 &
     & *(ZIPBL*PBUOM5(JL))**(ZCON2-1.0_JPRB)
  ENDIF
ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('VEXCSSAD_MOD:VEXCSSAD',1,ZHOOK_HANDLE)
END SUBROUTINE VEXCSSAD
END MODULE VEXCSSAD_MOD
